/*	///////////////////////////////////////////////////////////////////////
		Leaphy -- Likelihood estimation algorithms in phylogenetics
			Simon Whelan, University of Manchester

		Leaphy <data_file> <output_file> <Model> <InitialTree> <debug>

	Option definitions
	<data_file>		Any sequence data file in phylip or sequential format
	<output_file>	Output file
	<Model>			The model of sequence evolution (Name; +dG = gamma; +I = invariant sites
	<InitialTree>	The type of initial tree used for computation
					file=filename == tree from file
					all == best of starting trees; bionj = bionj; SA = stepwise addition
	<debug>			If debug is written then the program will run in debug mode (for checking routines are running properly)
	/////////////////////////////////////////////////////////////////////// */
#ifdef ISUNIX
#include "optimise.h"
#include "TreeList.h"
#include "model.h"
#include "interface.h"
#else
#include "../PhyloLib/optimise.h"
#include "../PhyloLib/TreeList.h"
#include "../PhyloLib/model.h"
#include "../PhyloLib/interface.h"
#endif

#if FUNC_COUNTERS == 1
	extern int Matrix_Log_Counter, MakeQ_Log_Counter, MakePT_Log_Counter, LFunc_Log_Counter, SubLFunc_Log_Counter, SPR_Log_Counter;
#endif
vector <double> GetPWVar(CBaseModel *Model, vector <double> *CurrentPW);
double GetDistVar(CBaseModel *M, CData *D,int Seq1, int Seq2, double CurrentDist);

/////////////////////////////////////////////////////////////////////////////////
// Main tree estimation routine
SBestModel DoTreeEstimation(CBaseModel *M, CData *D, CTree *T);

// Pointer to likelihood function that changes depending on the type of calculation performed
vector <STabuTree> TabuTrees;
vector <double> PWDists;
CPhyloDat PhyDat;
int TABU_RADIUS = DEFAULT_TABU_RADIUS, EXIT_OBS, OptObs;
double PROB_RAN_SEQ_REM = DEFAULT_PROB_RAN_SEQ_REM;
bool AllowPreOpt = true;

void DoInstructions();

bool WarningMulD;

// File specific global variables
int DebugOutput = 0;

int main(int argc, char *argv[])	{
	int i,NumModelReruns = 1;
	long RandomSeed = 0;
	string temp_string, Name;
	vector <string> Toks;
	TabuTrees.clear();
	WarningMulD = false;
	SBestModel Result;
	bool DoSurface = false, DoHardOptCheck = false;

	// Some initial verification
	if(min(LOOSE_RMSD_SUBSET,FULL_RMSD_SUBSET) < 105 && min(LOOSE_PARS_SUBSET,FULL_PARS_SUBSET) < 105) { Error("\nTrying to choose subsets of trees to examine based both on RMSD and parsimony"); }

	// Check what type of run it is
	if(argc > 5) {
		temp_string = argv[5];
		if(!strcmp(argv[5],"debug")) { DebugOutput = 1; }	// Run standard debug mode
		else if(temp_string.find("model") != string::npos) {
			PhyDat.NoCompress();
			Toks = Tokenise(temp_string,".");
			if(Toks.size() > 1) { NumModelReruns = atoi(Toks[1].c_str()); if(NumModelReruns < 0 || NumModelReruns > 100) {
				Error("\nNumber of model reruns is too large: " + int_to_string(NumModelReruns) + "\n\n"); }
			}
			if(NumModelReruns >1) {  DebugOutput = 2; } else { DebugOutput = 3; }
			PhyDat.SetModelAcc(FULL_LIK_ACC);
		} else if(temp_string.find("sim") != string::npos) { // Run simulation code
			Toks = Tokenise(temp_string,".");
			NumModelReruns = 1;
			if(Toks.size() > 1) { NumModelReruns = atoi(Toks[1].c_str()); if(NumModelReruns < 0) { Error("Trying to perform a negative number of simulations (" + int_to_string(NumModelReruns) + ")\n"); } }
			DebugOutput = 4;
			PhyDat.DoSimulation(argv[1], argv[2], NumModelReruns);
			cout << "\nDone!\n\n";
			exit(-1);
		}
		if(temp_string.find(".surface") != string::npos) { DoSurface = true; }
		if(temp_string.find(".optcheck") != string::npos || temp_string.find(".checkopt") != string::npos) { DoHardOptCheck = true; }
	}

	// Do data file
	if(argc > 1) {
		if(!strcmp(argv[1],"help")) { DoInstructions(); exit(-1); }
		PhyDat.SetIn(argv[1]); PhyDat.GetData();
	}
	// Do output file
	if(argc > 2) { PhyDat.SetOut(argv[2]); }
	if(argc > 3) {
		temp_string = argv[3];
		if(temp_string.find("file=") != string::npos) {
			temp_string.replace(0,5,"");
			PhyDat.ReadModel(temp_string);
			AllowPreOpt = false;
		//	PhyDat.LikelihoodSurfaces();
/*			cout.precision(12); cout << "\nHave read model: " << PhyDat.pModel()->lnL(true);
			cout << "\nModel again... " << PhyDat.pModel()->lnL(true);
			FOR(i,5) {
				cout << "\n["<<i<< "]: "<< PhyDat.pModel()->m_arL[i].LogP();

			}
//			cout << "\nPT"; PhyDat.pModel()->m_vpProc[0]->OutPT(cout,0);
	//		PhyDat.pModel()->ComplexOutput()	;
			cout << "\nModel:" << *PhyDat.pModel();
//			cout << "\nModel: " << PhyDat.pModel()->m_vpProc[0]->
			cout << "\nDone"; exit(-1);
	*/
			PhyDat.pModel()->ComplexOutput();
			cout << "\nModel:" << *PhyDat.pModel();
			exit(-1);
		} else {
			// Do model
			if(argc > 3) { PhyDat.SetModel(argv[3]); }
			// Do tree stuff
			if(argc > 4) { PhyDat.SetStartTree(argv[4]); }
			// Create the model and get it ready to use
			PhyDat.CreateModel();
	}	}

	///////////////////////////////////////////////////////
	// Get files and sequence data
	if(PhyDat.NeedMenu()) { ModelMenu(&PhyDat); PhyDat.CreateModel(); }
/*
	cout << "\nPrior to model set up: " << PhyDat.pModel()->lnL(true);
	cout << "\nDoing new data stuff..." << flush;
	CData *TEMP;
	TEMP = PhyDat.pData()->CreateAllPairs();
	cout << "\nDone" << flush;
	cout << *TEMP;
	exit(-1);
*/

//	cout << "\nDo bool AllowPreOpt = true;!!!"; system("sleep 3"); AllowPreOpt = false;


//	cout << "\nModel prep done" << flush;

//	cout.precision(12); system("sleep 1"); cout << "\nThe likelihood is: " << PhyDat.pModel()->lnL(); PhyDat.pTree()->OutBra(); PhyDat.pTree()->OutLabel(); cout << "\n" << *PhyDat.pTree(); exit(-1);

//	cout << "\nModel looks set up: " << PhyDat.pModel()->lnL(true) << flush;

#if DEVELOPER_BUILD == 1
		FOR(i,PhyDat.pTree()->NoBra()) {
			PhyDat.pTree()->SetStartCalc(i);
			cout << "\n\n>>>>>>>>>>>>>>>>>>> Doing lnL calc " << i << " <<<<<<<<<<<<<<<<<<<<<<<<<";
			cout.precision(12); cout << "\n\tLikelihood: " << PhyDat.pModel()->lnL();
		}
//		cout << "\nThe model is: " << *PhyDat.pModel()->m_vpProc[0];
//		exit(-1);

		cout << "\nThe P(t) matrices (WARNING: NEED TO RESET PERMISSIONS ON m_vpProc): ";
		FOR(i,3) {
			cout << "\n--- P("<< PhyDat.pTree()->B(i) << " --- ";
			int k,l;
			FOR(k,10) {
				cout << endl;
				FOR(l,10)	{
					cout << PhyDat.pModel()->m_vpProc[0]->PT(i)[(k*100)+l] << "  ";
			}	}
		}
#endif

/*
	// Some quick code that generates random trees
	cout << "\nDoing some random trees";
	FOR(i,1000000) {
		GetSAStartTree(PhyDat.pTree(), PhyDat.pModel(), false, randomtree, false);
		cout << "\nTree["<<i<<"]:\t" << *PhyDat.pTree();
	}
	exit(-1);
*/
	///////////////////////////////////////////////////////////////////////////////
	// Do general initialisation
	if(system("clear") == 1) { if(system("cls") == 1) { cout << "\nApologies for messy output -- unknown prompt..."; }; } ;
	cout.setf(ios::fixed); cout.precision(6);
	//////////////////////////////////////////////////////////////////////////////
	// Run type is determined by DebugOutput
	// ---
	// DebugOutput == 0: Standard tree estimation run
	// DebugOutput == 1: Debugging mode
	// DebugOutput == 2: Run a single model many times to attempt to avoid local optima
	// DebugOutput == 3: Run a single model once to calculate a likelihood
	// DebugOutput == 4: Do a simulation from a preloaded model
	// --
	// Set correct number of iterations
	if(DebugOutput == 2 || DebugOutput == 3) { if(DoHardOptCheck) { PhyDat.pModel()->OptNum(2000); } else { PhyDat.pModel()->OptNum(1000); } }
	// Do normal run
	if(DebugOutput == 0 || DebugOutput == 3) {		////////////////////////////////////////////////////////////////////////////////
		// Do the tree estimation
		////////////////////////////////////////////////////////////////////////////////
		Result = DoTreeEstimation(PhyDat.pModel(),PhyDat.pData(),PhyDat.pTree());

		////////////////////////////////////////////////////////////////////////////////
		// Output the exit details
		ofstream out(PhyDat.Out().c_str(), ios::app);
		out.setf(ios::fixed); out.precision(4);
		Result.m_vTree[0].OutBra(); Result.m_vTree[0].OutName();
		assert(Result.m_vTree.size() == 1);
		Result.m_vTree[0].OutLabel(); Result.m_vTree[0].OutName(); Result.m_vTree[0].OutBra();
		out << "\nBestTree\t" << Result.m_dlnL << "\t" << Result.m_vTree[0] << "\tObs=\t" << Result.m_iOptObs;
		out << "\n\n********************************************************************************\n*\tIf you use this tree in your research please cite:";
		out << "\n*\tS. Whelan. New Approaches to Phylogenetic Tree Search and their\n*\tApplication to Large Numbers of Protein Alignments.\n*\tSyst Biol. (2007) 56: 727-740";
		out << "\n********************************************************************************";
		out.close();
		PhyDat.pModel()->RestoreBestModel(&Result);
		if(PhyDat.pModel()->Tree()->Names().empty()) { PhyDat.pModel()->Tree()->SetNames(PhyDat.pModel()->m_pData->m_vsName); }
		cout << "\n\n~~~~~~ Exiting ~~~~~~\nBest likelihood=\t" << Result.m_dlnL << "\nBestTree\t" << *PhyDat.pModel()->Tree() << "\n";
		// Output the best model and its details to a seperate file
		string Name = PhyDat.Out() + "." + PhyDat.pModel()->Name() +".model";
		cout << "Outputting model to file <" << Name << ">";
		ofstream modout(Name.c_str());
		PhyDat.pModel()->ComplexOutput();
		if(PhyDat.Compress()) { PhyDat.pModel()->UncompressModel(); }
		modout << *PhyDat.pModel();
		modout.close();
		cout << " ... done";
		cout << "\n~~~~~~~~~~~~~~~~~~~~~\n";
	} else if(DebugOutput == 1) {	// Do basic debug run
	// Do debug run
		cout << "\nDoing model run...\n------------------------";
		if(PhyDat.pTree() == NULL) { cout << "\n\tCannot do debug likelihood when tree is provided..."; }

		if(AllowPreOpt) { cout << "\nDoing preoptimise..." << PreOpt(PhyDat.pModel(),false); }
		cout << "\n\nOriginal likelihood:   " << PhyDat.pModel()->lnL(true);
		cout << "\nOptimal likelihood:    " << flush; cout << FullOpt(PhyDat.pModel(),true,true,false) << flush << " == " << FullOpt(PhyDat.pModel(),true,true,false);
		// Check SA
		cout << "\nChecking SA tree:      " << flush;
		GetSAStartTree(PhyDat.pTree(),PhyDat.pModel(),false,lazy,DO_MP_SA);
		cout << FullOpt(PhyDat.pModel(),true,true) << flush;
		// Check bionj tree
		cout << "\nChecking bionj tree:   " << flush;
		PWDists = GetPW(PhyDat.pModel(),NULL,true);  // Get pairwise distances
		CTree T_bionj(DoBioNJ(PWDists,PhyDat.pData()->m_vsName,true),PhyDat.pData()->m_iNoSeq);
		*PhyDat.pModel()->m_pTree = T_bionj; PhyDat.pModel()->CleanFastCalc(true);
		cout << FullOpt(PhyDat.pModel(),true,true) << flush;


		cout << "\n\nSome debug stuff: ";
		PhyDat.pModel()->OutputDetail();
		cout << *PhyDat.pModel();
//		cout << "\nCalculations okay";
	} else if (DebugOutput == 2 && NumModelReruns > 1)	{	// Just run the model in extreme detail...
		CBaseModel *M = PhyDat.pModel();
		CTree OriTree; 	OriTree = *M->Tree();
		double lnL;
		int i,j,ModelNum=0;
		SBestModel BestModel; BestModel.m_dlnL = -BIG_NUMBER;
		cout << "\nDoing model only run... The model is: " << M->Name();
		if(!PhyDat.IsTree())	{ Error("\nCan only do model for a specified tree topology..."); }
		if(PhyDat.pModel()->Tree()->Names().empty()) { PhyDat.pModel()->Tree()->SetNames(PhyDat.pModel()->m_pData->m_vsName); }
		if(AllowPreOpt) { cout << "\nDoing preoptimise..." << PreOpt(PhyDat.pModel(),false); }

/*
PhyDat.CleanModel(); PhyDat.CleanTree();
PhyDat.ReadModel("Model.temp");
cout << "\nHave read model..." << flush;
M = PhyDat.pModel();
*/
		cout << "\nOriginal branches " << M->Tree()->Branches() << flush;
//		M->ComplexOutput();
//		cout << "\nInitial model: " << *M;
		lnL = M->lnL();
		cout << endl << M->Name() <<"_likelihood: " << lnL << flush;

		cout << "\nModel OptNum: " << M->OptNum();

		lnL = FullOpt(M,true,true,true,lnL,true,M->OptNum(),-BIG_NUMBER,PhyDat.ModelAcc(),false,true);
		cout << " -> " << lnL << flush;
		lnL = FullOpt(M,true,true,true,lnL,true,M->OptNum(),-BIG_NUMBER,PhyDat.ModelAcc(),false,true);
		cout << " -> " << lnL << flush;
		BestModel = M->ModelScores(lnL);
		cout << endl << M->Name() << "_par_names:\t" << M->ParNameOut();
		cout << endl << M->Name() << "_par_vals: \t" << M->ParValOut();
		cout << endl << M->Name() << "_Branches: \t" << M->Tree()->Branches();
		// Output the model details to a seperate file
		Name = PhyDat.Out() + "." + PhyDat.pModel()->Name() + "." + int_to_string(ModelNum++);
		cout << "\nOutputting to file: " << Name;
		ofstream mout(Name.c_str());
		M->ComplexOutput();
		mout << *M;
		mout.close();
		///////////////////////////////////////
		// Parameter shuffles
		cout << "\n\n------------------------------- Trying some shuffled parameters -----------------------------";
		cout << "\nRun#\tOrilnL\tOpt1\tOpt2\t" << M->ParNameOut(); FOR(j,M->Tree()->NoBra()) { cout << " Branch["<<j+1<<"]"; }
		ofstream Qout1("Q.output");
		Qout1 << "OriginalQ:" << flush; FOR(i,(int)M->m_vpProc.size()) { M->m_vpProc[i]->OutQ(-1,Qout1); }
		Qout1 << "\nPT: ";
		M->Tree()->SetB(0,0.25,true,true);
		M->m_vpProc[0]->Make_PT(0,true); M->m_vpProc[0]->OutPT(Qout1,0);
		Qout1 << "\nData: " << *M->m_pData;
		Qout1 << "\nDataEqm: " << M->m_pData->GetFreq();
		Qout1 << "\nEqm:	"; M->m_vpProc[0]->OutEqm(-1,Qout1);
		Qout1 << flush; Qout1.close();
//		cout << "\nThe Model: " << *M;
//		M->m_vpProc[0]->ShuffleParameters();
		FOR(i,NumModelReruns)	{
			M->RestoreBestModel(&BestModel);
			if(M->Name().find("THMM") == string::npos)	{ M->RandomiseParameters(true,true,true,true); }
			else { M->m_vpProc[0]->ShuffleParameters(); }
			*M->Tree() = OriTree; M->CleanFastCalc(true);
			cout << endl << i << "\t";
			lnL = M->lnL();
			cout  << lnL << "\t" << flush;
			lnL = FullOpt(M,true,true,true,lnL,true,M->OptNum(),-BIG_NUMBER,PhyDat.ModelAcc(),false,true);
			cout << lnL << "\t"<< flush;
			lnL = FullOpt(M,true,true,true,lnL,true,M->OptNum(),-BIG_NUMBER,PhyDat.ModelAcc(),false,true);
			cout << lnL << "\t"<< flush;
			cout << M->ParValOut() << " " << M->Tree()->Branches() << flush;
			// If the new model is the best, then store it
			if(lnL > BestModel.m_dlnL)	{
				BestModel = M->ModelScores(lnL);
			}
			// Output the model details to a seperate file
			Name = PhyDat.Out() + "." + M->Name() + "." + int_to_string(ModelNum++);
			ofstream mout2(Name.c_str());
			mout2 << "Likelihood:\t" << lnL << "\t";
			M->ComplexOutput();
			mout2 << *M;
			mout2.close();

//			ofstream Qout("Q.output",ios::app);
//			Qout << "\n----------------------------------------\nQ["<<i<<"]:"; FOR(i,(int)M->m_vpProc.size()) { M->m_vpProc[i]->OutQ(-1,Qout); } Qout.close();
//			cerr << "\nDone shuffle " << i;
		}
		// Output the best model and its details to a seperate file
		M->RestoreBestModel(&BestModel);
		// 1. Do the basics
		cout << endl << "best\t" << BestModel.m_dlnL << "\t" << BestModel.m_dlnL << "\t" << BestModel.m_dlnL << "\t";
		cout << M->ParValOut() << " " << M->Tree()->Branches() << flush;
		// 2. Do the file
		Name = PhyDat.Out() + "." + PhyDat.pModel()->Name() +".model";
		cout << "\nOutputting model to file <" << Name << ">";
		ofstream mout3(Name.c_str());
		M->ComplexOutput();
		mout3 << *M;
		mout3.close();
		cout << " done\n";
/*
		cout << "\n\tOriginal Likelihood: " << temp;
		// Do original optimise
		FullOpt(M,PhyDat.OptBra(),PhyDat.OptPar(),PhyDat.OptFre(),temp,true,DEFAULT_OPTNUM,-BIG_NUMBER,PhyDat.ModelAcc(),true);
		// Compress the model
		M->CompressModel();
		cout << "\n\tDone quick optimise and now compressing model" << flush;
		// Do a full optimise on the compressed model
		Best = M->ModelScores(FullOpt(M,PhyDat.OptBra(),PhyDat.OptPar(),PhyDat.OptFre(),temp,true,DEFAULT_OPTNUM,-BIG_NUMBER,FULL_LIK_ACC,true));
		cout << "\n\tOptimal likelihood for starting tree: " << Best.m_dlnL;
*/
	} else if (DebugOutput == 4) {
		cout << "\nAttempting to do a simulation should never reach here...";


		exit(-1);
	}
	if(DoHardOptCheck) {
		cout << "\nHard checking optimisation";
		Name = PhyDat.Out() + ".optcheck";
		ofstream optout(Name.c_str());
		DoHardOptCheck = CheckAllPar(PhyDat.pModel(),PhyDat.pModel()->lnL(true),PhyDat.pModel()->GetOptPar(),1.0E-3,optout);
		optout.close();
		cout << " -> <" << Name << "> "; if(DoHardOptCheck) { cout << " okay"; } else { cout << " fail"; }

	}
	if(DoSurface) {
		cout << "\nDoing likelihood surfaces";
		Name = PhyDat.Out() + ".surface";
		ofstream surout(Name.c_str());
		PhyDat.LikelihoodSurfaces(surout);
		surout.close();
		cout << " -> <" << Name << ">";
	}
#if FUNC_COUNTERS == 1
	cout << "\nFunction counters -- VMat(): " << Matrix_Log_Counter << ", lnL(): " << LFunc_Log_Counter << ", SublnL(): " << SubLFunc_Log_Counter << endl;
	OutputCounters(cout);
#endif
#if OUTPUT_TABULIST
	ofstream tabuout("Tabu.list");
	PhyDat.pModel()->DoParsimony();
	tabuout << TabuTrees.size() << "  " << PhyDat.NoSeq() << endl;
	tabuout << "#TreeNames\tTreeNums\tlnL\tPars\tRF2best\tIsOptima\n";
	FOR(i,(int)TabuTrees.size()) {
		TabuTrees[i].Tree.OutName(); TabuTrees[i].Tree.OutBra();
		tabuout<< TabuTrees[i].Tree << "\t";
		TabuTrees[i].Tree.NoOutName(); TabuTrees[i].Tree.NoOutBra();
		tabuout<< TabuTrees[i].Tree << "\t" << TabuTrees[i].lnL;
		// Do parsimony score
		PhyDat.pModel()->MakeNewTree(&TabuTrees[i].Tree,true);
		tabuout << "\t" << -PhyDat.pModel()->GetFullParsimony();
		// Do Robinson-Foulds distance to best
		tabuout << "\t" << Result.m_vTree[0].GetRFDist(TabuTrees[i].Tree);
		tabuout << "\t"; if(TabuTrees[i].Optima) { tabuout << "true"; } else { tabuout << "false"; }
		tabuout << endl;
	}
	tabuout.close();
#endif

	cout << "\n\nSuccessful exit\n";
	return 0;
}

/////////////////////////////////////////////////////////////////////////
// Tree estimation function
/////////////////////////////////////////////////////////////////////////
SBestModel DoTreeEstimation(CBaseModel *M, CData *D, CTree *T)	{
	int Rearrange, i,j;
	OptObs = 0; TABU_RADIUS = DEFAULT_TABU_RADIUS;
	bool SkipImpSamp = false, DoSPR, Output = true;
	double temp,FirstlnL;
	SBestModel Best, CurBest; Best.m_dlnL = -BIG_NUMBER;
	TabuTrees.clear();
	WarningMulD = false;
	string T_bionj;
	cout << "\nLeaphy: a phylogenetic tree searching program\n---------------------------------------------\n";
	cout << "\nInitialisation:";
	// Get tree lists
	cout << "\n\tGenerating lists of trees ... ";
	CTreeList FourSpTr(4); cout << " 4 ... " << flush;
	CTreeList FiveSpTr(5); cout << " 5 ... " << flush;
	CTreeList SixSpTr(6);  cout << " 6 done" << flush;

	// Set the correct likelihood accuracy for tree search
	if(DebugOutput != 3) { PhyDat.SetModelAcc(TREESEARCH_LIK_ACC); Output = false; }

	// Set resampling values
	int NUM_RESAMPLE = min(PhyDat.NoSeq(),50);
	EXIT_OBS = max(3 + (int) sqrt((double) NUM_RESAMPLE),3);
	if(EXIT_OBS > NUM_RESAMPLE) { EXIT_OBS = NUM_RESAMPLE; }
	/////////////////////////////////////////////////////
	// Deal with starting tree
	if(PhyDat.NoSeq() == 6 && DebugOutput != 3) {
		PhyDat.CleanTree();
		T = PhyDat.CreateTree("((1,2),(3,4),(5,6));");
		cout << "\n\tCreated 6 species tree for optimisation...";
		temp = M->lnL();
		// Do original optimise
		if(AllowPreOpt) { PreOpt(M,false); }
		FullOpt(M,PhyDat.OptBra(),PhyDat.OptPar(),PhyDat.OptFre(),temp,true,DEFAULT_OPTNUM,-BIG_NUMBER,PhyDat.ModelAcc(),true);
		// Compress the model
		if(PhyDat.Compress()) { M->CompressModel(); cout << "\n\tDone quick optimise and now compressing model" << flush; }
		// Do a full optimise on the compressed model
		Best = M->ModelScores(FullOpt(M,PhyDat.OptBra(),PhyDat.OptPar(),PhyDat.OptFre(),temp,true,DEFAULT_OPTNUM,-BIG_NUMBER,PhyDat.ModelAcc(),true));
		cout << "\n\tOptimal likelihood: " << Best.m_dlnL;
		// Get pairwise distances
		cout << "\n\tGetting pairwise distances ... " << flush;
		PWDists = GetPW(M,NULL,false); cout << " done" << flush;
	} else {
		if(!PhyDat.IsTree())	{
			//////////////////////// get a starting tree //////////////////////////////
			if(!PhyDat.SA() && !PhyDat.bionj()) { Error("No way defined to get starting tree...\n\n"); }
			T = PhyDat.CreateTree();
			if(PhyDat.SA()) {
				GetSAStartTree(T,M,true,lazy,DO_MP_SA);
				if(DO_MP_HILLCLIMB == 1 && DO_MP_SA == 1) { DoMPHillClimb(M,&FourSpTr,&FiveSpTr,&SixSpTr); }
				temp = M->lnL();
				cout << "\n\tOriginal Likelihood: " << temp << " ... doing initial optimise (this may take some time)" <<flush;
				// Do original optimise
				if(AllowPreOpt) { PreOpt(M,false); }
				FullOpt(M,PhyDat.OptBra(),PhyDat.OptPar(),PhyDat.OptFre(),temp,true,DEFAULT_OPTNUM,-BIG_NUMBER,PhyDat.ModelAcc(),Output);
				// Compress the model
				if(PhyDat.Compress()) { M->CompressModel();	cout << "\n\tDone quick optimise and now compressing model" << flush; }
				// Do a full optimise on the compressed model
				Best = M->ModelScores(FullOpt(M,PhyDat.OptBra(),PhyDat.OptPar(),PhyDat.OptFre(),temp,true,DEFAULT_OPTNUM,-BIG_NUMBER,PhyDat.ModelAcc(),Output));
				FirstlnL = Best.m_dlnL;
				cout << "\n\tOptimal likelihood for stepwise addition tree: " << Best.m_dlnL;
				// Get pairwise distances
				cout << "\n\tGetting pairwise distances ... " << flush;
				PWDists = GetPW(M,NULL,false); cout << " done" << flush;
			} else { // Need to get pairwise distances if no stepwise addition
				cout << "\n\tGetting pairwise distances ... " << flush;
				PWDists = GetPW(M,NULL,true); cout << " done" << flush;
			}
			if(PhyDat.bionj())	{
				cout << "\n\n------------------------- Preparing bionj tree ---------------------------------\n";
				// Get the bionj tree
				CTree T_bionj(DoBioNJ(PWDists,D->m_vsName,true),D->m_iNoSeq,false,D);
				M->MakeNewTree(&T_bionj,true);
				M->CleanFastCalc(true);
				if(DO_MP_HILLCLIMB == 1) { DoMPHillClimb(M,&FourSpTr,&FiveSpTr,&SixSpTr); }
				temp = M->lnL();
				cout << "\n\tOriginal Likelihood: " << temp << flush;
				if(PhyDat.SA()) { // If already done some optimisation then do quicky
					temp = FullOpt(M,PhyDat.OptBra(),PhyDat.OptPar(),PhyDat.OptFre(),temp,true,DEFAULT_OPTNUM,-BIG_NUMBER,PhyDat.ModelAcc(),Output); cout << "\n\tbionj tree optimal likelihood: " << temp;
					if(temp > Best.m_dlnL) {
						Best = M->ModelScores(temp); cout << " ... better than SA" << flush;
					} else { cout << " ... reverting to SA" << flush; M->RestoreBestModel(&Best); }
				} else { // Do the quick model optimise and reduce it
					// Do original optimise
					if(AllowPreOpt) { PreOpt(M,false); }
					FullOpt(M,PhyDat.OptBra(),PhyDat.OptPar(),PhyDat.OptFre(),temp,true,DEFAULT_OPTNUM,-BIG_NUMBER,PhyDat.ModelAcc(),Output);
					// Compress the model
					if(PhyDat.Compress()) { M->CompressModel(); cout << "\n\tDone quick optimise and now compressing model" << flush; }
					// Do a full optimise on the compressed model
					Best = M->ModelScores(FullOpt(M,PhyDat.OptBra(),PhyDat.OptPar(),PhyDat.OptFre(),temp,true,DEFAULT_OPTNUM,-BIG_NUMBER,PhyDat.ModelAcc(),Output));
					FirstlnL = Best.m_dlnL;
					cout << "\n\tOptimal likelihood for bionj tree: " << Best.m_dlnL;
					cout << "\n\tUpdating pairwise distances ... " << flush;
					PWDists = GetPW(M,&PWDists,false); cout << " done" << flush;
		}	}	} else {
			//////////////////////// Do an initial tree //////////////////////////////
			cout << "\n\n------------------------- Preparing original tree ------------------------------\n" << flush;
			cout << "\nCall lnL" << flush;
			temp = M->lnL();
			cout << " ... done";
			cout << "\nDone likelihood " << temp << flush;
			// Do original optimise
			if(AllowPreOpt) {
				cout << "\nPerforming preoptimisation (may take some time) from " << temp << flush;
				PreOpt(M,false);
				cout << " done" << flush;
				temp = FullOpt(M,PhyDat.OptBra(),PhyDat.OptPar(),PhyDat.OptFre(),temp,true,DEFAULT_OPTNUM,-BIG_NUMBER,PhyDat.ModelAcc(),Output);
				// Compress the model
				if(PhyDat.Compress()) { M->CompressModel(); cout << "\n\tDone quick optimise and now compressing model" << flush; }
			}
			cout << "\nPerforming full optimisation from " << temp;
			// Do a full optimise on the compressed model
			Best = M->ModelScores(FullOpt(M,PhyDat.OptBra(),PhyDat.OptPar(),PhyDat.OptFre(),temp,true,M->OptNum(),-BIG_NUMBER,PhyDat.ModelAcc(),Output));
			FirstlnL = Best.m_dlnL;
			cout << "\n\tOptimal likelihood for starting tree: " << Best.m_dlnL;
			// Get pairwise distances
			if(DebugOutput != 3)	{
				cout << "\n\tGetting pairwise distances ... " << flush;
				PWDists = GetPW(M,NULL,false); cout << " done" << flush;
	}	}	}
	cout << "\n\n----------------------------------- Done ---------------------------------------\n";
	// Set as first tabu
	if(IsTabu(T,false,true,M,Best.m_dlnL) == true) { cout << "\nImpossible: tree thought tabu at start..."; exit(-1); }
/*
////////////////////////// Some ratchetting stuff //////////////////////////////
const int NumIter = 100;
int i,count = 0; T->NoOutBra();
CTree *TempT, TStore; TStore = *M->Tree();
// Do ratchet bionj trees
ofstream r_bionj("ratchet_bionj");
r_bionj << "Tree\tOri_lnL\tOpt_lnL\n";
cout << "\n----- Doing ratchets of bionj -----\n";
FOR(i,NumIter)	{
	if(i%50==0) { cout << "\n" << flush; } cout << "." << flush;
	TempT = new CTree(DoRandomBioNJ(M,D),M->NoSeq());
	*M->Tree() = *TempT;
	r_bionj << *TempT << "\t" << M->lnL() << "\t" << flush;
	r_bionj << FullOpt(M) << endl << flush;
	*M->Tree() = TStore;
	delete TempT; TempT = NULL;
}
r_bionj.close();
// Do fast SA trees
D->CleanBootstrap();		// Just to make sure
ofstream SAout("SA_normal");
SAout << "Tree\tOri_lnL\tOpt_lnL\n";
cout << "\n\n----- Doing SA normal -----\n";
TempT = new CTree;
FOR(i,NumIter)	{
	if(i%50==0) { cout << "\n" << flush; } cout << ".";
	GetSAStartTree(M->Tree(),M,false,true,DO_ML_SA);
	M->CleanFastCalc(true);
	SAout << *M->Tree() << "\t" << M->lnL() << "\t" << flush;
	SAout << FullOpt(M) << endl << flush;
}
SAout.close();
// Do ratchet fast SA trees
// Do fast SA trees
D->CleanBootstrap();		// Just to make sure
ofstream ratchetSA("ratchet_SA");
ratchetSA<< "Tree\tOri_lnL\tOpt_lnL\n";
cout << "\n\n----- Doing SA ratchet -----\n";
FOR(i,NumIter)	{
	if(i%50==0) { cout << "\n" << flush; } cout << ".";
	D->Bootstrap();
	GetSAStartTree(M->Tree(),M,false,true,DO_ML_SA);
	D->CleanBootstrap();
	M->CleanFastCalc(true);
	ratchetSA<< *M->Tree() << "\t" << M->lnL() << "\t" << flush;
	ratchetSA<< FullOpt(M) << endl << flush;
}
delete TempT; TempT = NULL;
ratchetSA.close();
exit(-1);
//////////////////////////////////////////////////////////////////////////
*/
	////////////////////////////////////////////////////////////////////////////////////////
	// Stop everything here if not doing tree search
	if(DebugOutput == 3) {
		// 2. Do the file
		string Name;
		Name = "Model." + M->Name() + ".best";
		ofstream mout3(Name.c_str());
		M->ComplexOutput();
		mout3 << *M;
		mout3.close();
		return Best;
	}

	// Ensure tree estimation doesn't happen for parameters deccribing pre-specified groups of branches
	if(PhyDat.pTree()->BranchLabels()) { Error("\nTrying to estimate a tree with parameters varying in pre-defined groups\n"); }

	cout << "\nTo here... 2" << flush;

	///////////////////////////////////////////////////////////////////////////////////////////
	// Output the pairwise distances...
	/////////////////////////////////////////////////////////////////////////////////////////
	string pairname = PhyDat.Out();
	if(pairname.find(".output") != string::npos) {
		pairname = pairname.replace(pairname.find(".output"),7,"");
	}
	pairname += ".dist";
	ofstream distout(pairname.c_str());
	distout << PhyDat.NoSeq();
	double MaxDist = 0.0, MinDist = BIG_NUMBER;
	FOR(i,PhyDat.NoSeq()) {
		distout << endl;
		FOR(j,PhyDat.NoSeq())	{
			distout << PWDists[(i*PhyDat.NoSeq())+j] << "\t";
			if(i==j) { continue; }
			if(PWDists[(i*PhyDat.NoSeq())+j] > MaxDist) { MaxDist = PWDists[(i*PhyDat.NoSeq())+j]; }
			if(PWDists[(i*PhyDat.NoSeq())+j] < MinDist) { MinDist = PWDists[(i*PhyDat.NoSeq())+j]; }
		}
	}
	distout << "\nMaximum_distan-- Donece: " << MaxDist << "\nMinimum_distance: " << MinDist;
	distout.close();

	///////////////////////////////////////////////////////////////////////////////////////////
	// This samples tree space and hill-climbs
	///////////////////////////////////////////////////////////////////////////////////////////
	// vectors holding the optimal likelihoods and the number of times they've been seen
	// Will also try larger rearrangements
	ofstream initout(PhyDat.Out().c_str());
	initout << "StartingTree:\t" << *T << "\t=\t" << Best.m_dlnL;
	initout.close();


	FOR(Rearrange,NUM_RESAMPLE)	{
		cout << "\n\n\t=======================================================================================";
		cout << "\n\t      Rearrangement " << Rearrange+1 << " / " << NUM_RESAMPLE << " (T_RAD = " << TABU_RADIUS * 2;
		if(OptObs >= EXIT_OBS - COMP_SNAP_AWAY_2_OPT_OBS) { cout << ":ForceAway"; }
		cout << ")\tOptObs " << OptObs+1 << " / " << EXIT_OBS << "\tExamined " << TabuTrees.size() << " tree(s)";

		cout << "\n\t=======================================================================================";
		cout << "\n\nBest likelihood: " << Best.m_dlnL;
		cout << "\nBest Tree:     " << Best.m_vTree;
		cout << "\nProposed Tree: " << *M->m_pTree << flush;


		ofstream out(PhyDat.Out().c_str(), ios::app); out << "\nSWitch iteration " << Rearrange; out.close();
		if(Rearrange > 0 && IsTabu(M->m_pTree,true,true,M,M->lnL()) == true && DO_SNAP == 1) {
			// Catch for when the optimiser screws up
			if(WarningMulD == true) {
				WarningMulD = false;
				if(CurBest.m_dlnL > Best.m_dlnL + FULL_LIK_ACC)	{ Best = CurBest; }
			} else { cout << "\nProposed tree tabu..."; exit(-1); }
		}
#if DO_SNAP == 1
		if(DO_SPR_FULL == 1) { Error("\nCannot do SNAP and full SPR search..."); }
		////////////////////////////////////////////
		// Perform SNAP optimisation
		////////////////////////////////////////////
		CurBest = M->ModelScores(TreeSNAP(&FourSpTr,&FiveSpTr,&SixSpTr,M,true,lazy));

		cout << "\nBest likelihood after SNAP: " << CurBest.m_dlnL;

		DoSPR = false;
		if(CurBest.m_dlnL < 0)	{	// If this tree hasn't been seen before
			// Store if this beats what its seen previously
			if(CurBest.m_dlnL > Best.m_dlnL + FULL_LIK_ACC)	{
				Best = CurBest;
				// Reset OptObs only if it is a different tree topology...
				if(MinTabuRFDist(&CurBest.m_vTree[0],true) > 0) { TABU_RADIUS = DEFAULT_TABU_RADIUS; PROB_RAN_SEQ_REM = DEFAULT_PROB_RAN_SEQ_REM; OptObs = 0; DoSPR = true; }

			}
			else { M->RestoreBestModel(&Best); }
		}
		// Output current best
		ofstream fout(PhyDat.Out().c_str(), ios::app);
		assert(Best.m_vTree.size() == 1);
		Best.m_vTree[0].OutBra();
		fout.precision(6); fout.setf(ios::fixed);
		Best.m_vTree[0].OutName();
		fout << "\nCurOpt:\t" << Best.m_vTree[0] << "\t" << Best.m_dlnL << "\tOptObs=\t" << OptObs;
		Best.m_vTree[0].NoOutName();
		Best.m_vTree[0].NoOutBra();
		fout.close();

		////////////////////////////////////////////////////////
		// Perform a round of SPR if a new tree is discovered
		////////////////////////////////////////////////////////
#if DO_SPR_SECOND == 1
		double NewlnL;
		if(DoSPR) {
			M->RestoreBestModel(&Best);
			cout << "\n\n------------------ Trying SPR: saves sticking around looking at bad trees ------------------";
			cout << "\nInitial tree: " << *M->Tree();
			NewlnL = TreeSPR(M,true,Best.m_dlnL,true,true,10,lazy);
			ofstream sprout(PhyDat.Out().c_str(), ios::app);
			if(NewlnL > Best.m_dlnL + FULL_LIK_ACC) {

				cout << "\nBest: " << Best.m_vTree[0];
				cout << "\nCurr: " << *M->Tree();
				cout << "\nMinTabuRFDist: " << MinTabuRFDist(&Best.m_vTree[0],false);

				if(!IsSameTree(&Best.m_vTree[0],M->Tree()))	{
					cout << "\nSPR found improved tree! lnL = " << NewlnL;
					Best = M->ModelScores(NewlnL); OptObs = 0; TABU_RADIUS = DEFAULT_TABU_RADIUS; PROB_RAN_SEQ_REM = DEFAULT_PROB_RAN_SEQ_REM;
					SkipImpSamp = true;
					sprout << "\nSPR gives better tree:\t" << *M->Tree() << "\t" << CurBest.m_dlnL;
					OptObs--;
					Rearrange--;
				} else { cout << "\nSPR improved likelihood..."; Best = M->ModelScores(FullOpt(M)); }
			} else {
				sprout << "\nSPR gives no improvement";
			}
			sprout.close();
		}
		M->RestoreBestModel(&Best);
#endif
#elif DO_SPR_FULL == 1

		CurBest = M->ModelScores(TreeSPR(M,true,M->lnL(),true,true,100,lazy));
		// Do the first tree if required
//		if(Rearrange <= 0 && Best.m_dlnL - FirstlnL > 10.0 && M->Compressed()) {
		if(Rearrange <= 0) {
			ofstream upout(PhyDat.Out().c_str(), ios::app);
			M->RestoreBestModel(&CurBest);
			cout << "\nUncompressing original model and reoptmising so likelihoods are more accurate... " << Best.m_dlnL << flush;
			upout << "\nReoptimised model and tree parameters: " << Best.m_dlnL;
			M->UncompressModel();
			Best = M->ModelScores(FullOpt(M,PhyDat.OptBra(),PhyDat.OptPar(),PhyDat.OptFre(),-BIG_NUMBER,true,DEFAULT_OPTNUM,-BIG_NUMBER,PhyDat.ModelAcc(),false));
			TabuTrees.clear();
			cout << " -> " << Best.m_dlnL << flush;
 			upout << " -> " << Best.m_dlnL << flush;
			M->LockModel();
			cout << " ... locked";
			// Do branch optimisation to high degree of accuracy
			cout << "\nBestlnL: " << Best.m_dlnL;
			M->GetOptPar(true,true,false,false);
			Best = M->ModelScores(M->FastBranchOpt(-BIG_NUMBER,FULL_LIK_ACC));
			cout << " --> " << Best.m_dlnL;
			upout.close();
			// Set as first tabu
			if(IsTabu(T,false,true,M,Best.m_dlnL) == true) { cout << "\nImpossible: tree thought tabu at start..."; exit(-1); }
			cout << "\n\tUpdating pairwise distances ... " << flush;
			PWDists = GetPW(M,&PWDists,false); cout << " done" << flush;
		} else {
			// Otherwise normal tabu search
			// Optimise branches to high degree of accuracy
			M->GetOptPar(true,true,false,false);
			CurBest = M->ModelScores(M->FastBranchOpt(-BIG_NUMBER,FULL_LIK_ACC));
			if(CurBest.m_dlnL < 0)	{	// If this tree hasn't been seen before
				// Store if this beats what its seen previously
				if(CurBest.m_dlnL > Best.m_dlnL + FULL_LIK_ACC)	{
					Best = CurBest;
					// Reset OptObs only if it is a different tree topology...
					if(MinTabuRFDist(&CurBest.m_vTree[0],true) > 0) { TABU_RADIUS = DEFAULT_TABU_RADIUS; PROB_RAN_SEQ_REM = DEFAULT_PROB_RAN_SEQ_REM; OptObs = 0; DoSPR = true; }
				}
				else { M->RestoreBestModel(&Best); }
			}
		}
		// Output current best
		ofstream fout(PhyDat.Out().c_str(), ios::app);
		assert(Best.m_vTree.size() == 1);
		Best.m_vTree[0].OutBra();
		fout.precision(6); fout.setf(ios::fixed);
		Best.m_vTree[0].OutName();
		fout << "\nCurOpt:\t" << Best.m_vTree[0] << "\t" << Best.m_dlnL << "\tOptObs=\t" << OptObs;
		Best.m_vTree[0].NoOutName();
		Best.m_vTree[0].NoOutBra();
		fout.close();
#endif

/////////////////////////////////////////////////////

		/////////////////////////////////////////////
		// Exit criteria
		/////////////////////////////////////////////

		OptObs++;
		if(OptObs > 3) { TABU_RADIUS = min(TABU_RADIUS+2,(int)min((double)PhyDat.NoSeq() * 0.5,25)); }
		if(PhyDat.NoSeq() <= 6 || OptObs == EXIT_OBS || NUM_RESAMPLE == 1) { break; }
#if DO_IMPSAMPLE == 1
		/////////////////////////////////////////////
		// Resampling of tree space
		// Define number of attempts dependent on OptObs
		if(Rearrange < NUM_RESAMPLE-1 && !SkipImpSamp) {
			AgeTabuTrees();
			CurBest = M->ModelScores(GetTree(M,&FourSpTr,&FiveSpTr,&SixSpTr));
		}
		SkipImpSamp = false;
#else
		break;
#endif
#if DO_SNAP == 0
		if(CurBest.m_dlnL > Best.m_dlnL + FULL_LIK_ACC)	{
			Best = CurBest;
			OptObs = 0;
		}
#endif
	}
	M->RestoreBestModel(&Best);
	Best.m_iOptObs = OptObs;
	return Best;
}

// A function that explains how to run the program from the command line
void DoInstructions()   {
        cout << "\nTo run Leaphy from the commandline:\n\n\t./Leaphy  <data_file> <output_file> <Model> <Start_tree> <Advanced_options>\n";
        cout << "\n\t<data_file>:  \tThe input file in phylip or sequential format ([advanced] or the model file if simlulating data)";
        cout << "\n\t<output_file>:\tThe location of the output file ([advanced] sequence alignments if simulating data)";
        cout << "\n\t<Model>:      \tThe model";
        cout << "\n\t<Start_tree>: \tWhere the tree estimation program starts\n\t\t\"all\": Try stepwise addition and bionj";
        cout << "\n\t\t\"bionj\": use only bionj\n\t\t\"SA\": use only stepwise addition";
        cout << "\n\t\t\"file=XXXX\": Load an input tree from file XXXX ([advanced] will be ignored for simulated data)";
        cout << "\n\t<Advanced_options>:\n\t\t\tmodel = run only the model and do not optimise the tree";
        cout << "\n\t\t\tmodel.n = run only the model, but do so n times in an attempt to avoid local optima";
        cout << "\n\t\t\tsim.n = simulate n sets of sequence alignments. If .n isn't included Leaphy will simulate a single alignment";
        cout << "\n\n";
}
